Full Testate Caribou and Hole Bog Reconstructions
John Tipton
6/8/2017
set.seed(11)
library(BayesComposition)
library(knitr)
library(ggplot2)
library(gjam)
library(rioja)
library(analogue)
N <- 500
N_pred <- 25
d <- 4
n_knots <- 30Load Testate Data and R code - Full Data
testatePlotData <- data.frame(species=as.factor(rep(colnames(y), each=N)),
Count=c(as.matrix(y_prop)),
Wetness=rep(X, times=dim(y)[2]))
ggplot(testatePlotData, aes(x=Wetness, y=Count, color=species, group=species)) +
geom_point(alpha=0.25) +
# geom_line(stat="smooth", method="loess", aes(y=Count, x=Wetness),
# alpha=0.5, lwd=1.25) +
theme(legend.position="none") + ggtitle("Testate Composition vs. Water Table Depth") +
labs(x="Water Table Depth", y="Composition") +
theme(plot.title=element_text(size=24, face="bold", hjust=0.5)) +
theme(axis.text.x = element_text(size = 22),
axis.text.y = element_text(size = 22),
axis.title.x = element_text(size = 22),
axis.title.y = element_text(size = 22))
y <- as.matrix(y)
mean_X <- mean(X)
sd_X <- sd(X)
# X <- (X - mean_X) / sd_X
X_knots <- seq(min(X, na.rm=TRUE)-1.25*sd(X, na.rm=TRUE),
max(X, na.rm=TRUE)+1.25*sd(X, na.rm=TRUE), length=n_knots)
params <- list(n_adapt=10000, n_mcmc=20000, n_thin=20,
likelihood="dirichlet-multinomial",
function_type = "basis",
additive_correlation=FALSE,
multiplicative_correlation=FALSE,
message=100, n_chains=4, n_cores=4, df=8)
save_directory <- "./mvgp-test/"
save_file <- "mvgp-dm-testate-bspline.RData"
progress_directory <- "./mvgp-test/"
progress_file <- "mvgp-dm-testate-bspline.txt"
if (file.exists(paste0(save_directory, save_file))) {
## check if MCMC output exists
load(paste0(save_directory, save_file))
} else {
## potentially long running MCMC code
out <- fit_compositional_data(y=y, X=X, params=params,
progress_directory = "./mvgp-test/",
progress_file = "mvgp-dm-testate-bspline.txt",
save_directory = save_directory,
save_file = save_file)
}## Loading required package: snow## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('/Users/Tips/BayesComposition/fit-mvgp-dirichlet-
## multinomial-testate.Rmd',~+~~+~encoding~+~## R Version: R version 3.4.1 (2017-06-30)## snowfall 1.84-6.1 initialized (using snow 0.4-2): parallel execution on 4 CPUs.## Loading required namespace: rlecuyer## Library BayesComposition loaded.## Library BayesComposition loaded in cluster.## Library coda loaded.## Library coda loaded in cluster.##
## Stopping cluster## extract posterior samples
samples <- extract_compositional_samples(out)
alpha_post <- samples$alpha
beta_post <- samples$betaRhat <- make_gelman_rubin(out)
layout(matrix(1:3, 3, 1))
# hist(Rhat[grepl("alpha=", names(Rhat))], main = "Rhat for alpha")
hist(Rhat[grepl("beta", names(Rhat))], main = "Rhat for beta")
hist(Rhat, main="All parameters")
Rhat[!is.finite(Rhat)] <- NA
max(unlist(na.omit(Rhat)))## [1] 1.026868
alpha_post_mean <- apply(alpha_post, c(2, 3), mean)
## force the sum to one constraint
p_alpha_post_mean <- alpha_post_mean
for (i in 1:N) {
p_alpha_post_mean[i, ] <- p_alpha_post_mean[i, ] / sum(p_alpha_post_mean[i, ])
}
y_percentages <- y
for (i in 1:N) {
y_percentages[i, ] <- y_percentages[i, ] / sum(y_percentages[i, ])
}
fitPlotData <- data.frame(species=as.factor(rep(1:d, each=N)),
count=c(y_percentages),
depth=rep(X, times=d),
alpha=c(p_alpha_post_mean))
g1_post <- ggplot(fitPlotData, aes(x=depth, y=count, color=species, group=species)) +
geom_point(alpha=0.25) + theme(legend.position="none") +
ggtitle("MVGP vs. depth") +
geom_line(aes(x=depth, y=alpha, col = species), fitPlotData, lwd=1.25)
g2_post <- ggplot(fitPlotData, aes(x=depth, y=count, color=species, group=species)) +
geom_point(alpha=0.25) + theme(legend.position="none") +
ggtitle("MVGP vs. depth by species") +
geom_line(aes(x=depth, y=alpha, col = species), fitPlotData, lwd=1.25) +
facet_wrap( ~ species, ncol = 4)
multiplot(g1_post, g2_post, cols=2)## Loading required package: grid
load core for reconstruction
##
## Reconstruction period data at Hole Bog
##
fileToLoad <- '~/testate/data/Hole Bog 2005 Core, raw data.csv'
n = max(count.fields(fileToLoad, sep=","), na.rm=TRUE)
xx = readLines(fileToLoad)
xx = xx[-c(1:6, 47:56)]
.splitvar = function(x, sep, n) {
var = unlist(strsplit(x, split=","))
length(var) = n
return(var)
}
xx = do.call(cbind, lapply(xx, .splitvar, sep=",", n=n))
xx = apply(xx, 1, paste, collapse=",")
y_recon = read.csv(text=xx, sep=",", header=TRUE)
## processing data to have names first 3 letters of genus and species
genus_hole <- tolower(substr(names(y_recon), 1, 3))
species_hole <- substr(sub(".*\\.", "", sub("\\.catinus", "", sub("\\.minor", "", sub("\\.\\.", "", sub("\\.type", "", names(y_recon)))))), 1, 3)
names(y_recon) <- paste(genus_hole, species_hole, sep="")
source("~/testate/data/join-testate-caribou.R")
##
## align names between data sets
## get expert input here!!
## We are removing many species
##
N_recon_full <- sum(y_recon)
idx1 <- colnames(y) %in% names(y_recon)
# names(y_recon)[!idx1]
# names(y_cal)[!idx2]
## add in any species in calibation data not in reconstruction data
for (i in colnames(y)[!idx1]){
y_recon[, i] <- rep(0, dim(y_recon)[1])
}
idx2 <- names(y_recon) %in% colnames(y)
## remove species in reconstruction data not in calibration data
for (i in names(y_recon)[!idx2]){
y_recon[, i] <- NULL
}
y_recon <- y_recon[, order(colnames(y_recon))]
## check the names and orderings
all.equal(colnames(y), colnames(y_recon))## [1] TRUE## Note: we are reducing the testate counts by about 25%...
N_recon_reduced <- sum(y_recon)
N_recon_full - N_recon_reduced## [1] 903pred <- predict_compositional_data(as.matrix(y_recon), X,
params=params, samples=samples,
progress_directory=progress_directory,
progress_file = "DM-predict.txt")## Starting predictions, running for 4000 iterations
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Prediction iteration 100
## Prediction iteration 200
## Prediction iteration 300
## Prediction iteration 400
## Prediction iteration 500
## Prediction iteration 600
## Prediction iteration 700
## Prediction iteration 800
## Prediction iteration 900
## Prediction iteration 1000
## Prediction iteration 1100
## Prediction iteration 1200
## Prediction iteration 1300
## Prediction iteration 1400
## Prediction iteration 1500
## Prediction iteration 1600
## Prediction iteration 1700
## Prediction iteration 1800
## Prediction iteration 1900
## Prediction iteration 2000
## Prediction iteration 2100
## Prediction iteration 2200
## Prediction iteration 2300
## Prediction iteration 2400
## Prediction iteration 2500
## Prediction iteration 2600
## Prediction iteration 2700
## Prediction iteration 2800
## Prediction iteration 2900
## Prediction iteration 3000
## Prediction iteration 3100
## Prediction iteration 3200
## Prediction iteration 3300
## Prediction iteration 3400
## Prediction iteration 3500
## Prediction iteration 3600
## Prediction iteration 3700
## Prediction iteration 3800
## Prediction iteration 3900
## Prediction iteration 4000n_samples <- length(pred$X[, 1])
N_pred <- length(pred$X[1, ])
X_ci <- apply(pred$X, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ]
hole.df <- data.frame(
Covariate=c(X_ci),
Observation=factor(rep((1:N_pred), each=n_samples*0.95))
)
ggplot(hole.df, aes(Observation, Covariate)) +
geom_violin(position="identity") +
scale_x_discrete(breaks=seq(5, 25, 5)) +
labs(x="Observation", y="Unobserved July Temperature")
n_samples <- length(pred$X[, 1])
N_pred <- length(pred$X[1, ])
preds_WA <- matrix(0, n_samples, N_pred)
preds_MLRC <- matrix(0, n_samples, N_pred)
preds_MAT <- matrix(0, n_samples, N_pred)
# preds_GJAM <- matrix(0, n_samples, d)
for (i in 1:N_pred) {
preds_WA[, i] <- rnorm(n_samples, pred_mu_WA[i], pred_sd_WA[i])
preds_MLRC[, i] <- rnorm(n_samples, pred_mu_MLRC[i], pred_sd_MLRC[i])
preds_MAT[, i] <- rnorm(n_samples, pred_mu_MAT[i], pred_sd_MAT[i])
# preds_GJAM[, j] <- rnorm(n_samples, pred_mu_GJAM, pred_sd_GJAM)
}
X_ci <- rbind(apply(pred$X, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_WA, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_MLRC, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_MAT, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_rf, 1, sort)[(0.025*n_samples+1):(0.975*n_samples), ])
# apply(preds_GJAM, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ])
hole.df <- data.frame(
Covariate=c(X_ci),
Observation=factor(rep((1:N_pred), each=5*n_samples*0.95)),
Model=rep(c("MVGP", "WA", "MLRC", "MAT", "RF"), each=n_samples*0.95)
)
ggplot(hole.df, aes(Observation, Covariate, color=Model)) +
geom_violin(position="dodge") +
scale_x_discrete(breaks=seq(5, 25, 5)) +
labs(x="Observation", y="Unobserved Water Depth") +
facet_wrap(~ Model, ncol=2) +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model))
ggplot(hole.df, aes(Observation, Covariate, color=Model)) +
scale_x_discrete(breaks=seq(5, 25, 5)) +
labs(x="Observation", y="Unobserved Water Depth") +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
stat_summary(fun.ymin = function(z) {quantile(z, 0.025)},
fun.ymax = function(z) {quantile(z, 0.975)},
geom = "ribbon", aes(Observation, Covariate,
group=Model, color=Model, fill=Model), alpha=0.25) +
ylim(-10, 60)## Warning: Removed 2148 rows containing non-finite values (stat_summary).
## Warning: Removed 2148 rows containing non-finite values (stat_summary).
ggplot(hole.df, aes(Observation, Covariate, color=Model)) +
scale_x_discrete(breaks=seq(5, 25, 5)) +
labs(x="Observation", y="Unobserved Water Depth") +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
stat_summary(fun.ymin = function(z) {quantile(z, 0.025)},
fun.ymax = function(z) {quantile(z, 0.975)},
geom = "ribbon", aes(Observation, Covariate,
group=Model, color=Model, fill=Model), alpha=0.25) +
facet_wrap(~Model, ncol=2) +
ylim(-10, 60)## Warning: Removed 2148 rows containing non-finite values (stat_summary).
## Warning: Removed 2148 rows containing non-finite values (stat_summary).
Caribou bog
# fileToLoad <- read.csv("~/Documents/connoR/Caribou/cariboucounts.csv")
fileToLoad <- read.csv("~/testate/data/rawData/cariboucounts.csv")
y_recon <- fileToLoad[,10:60]
#source("~/Documents/connoR/testate/data/join-testate-caribou.R")
N_recon_full <- sum(y_recon)
idx1 <- colnames(y) %in% names(y_recon)
# names(y_recon)[!idx1]
# names(y_cal)[!idx2]
## add in any species in calibation data not in reconstruction data
for (i in colnames(y)[!idx1]){
y_recon[, i] <- rep(0, dim(y_recon)[1])
}
idx2 <- names(y_recon) %in% colnames(y)
## remove species in reconstruction data not in calibration data
for (i in names(y_recon)[!idx2]){
y_recon[, i] <- NULL
}
y_recon <- y_recon[, order(colnames(y_recon))]
## check the names and orderings
all.equal(colnames(y), colnames(y_recon))## [1] TRUE## Note: we are reducing the testate counts by about 25%...
N_recon_reduced <- sum(y_recon)
N_recon_full - N_recon_reduced## [1] 9046pred <- predict_compositional_data(as.matrix(y_recon), X,
params=params, samples=samples,
progress_directory=progress_directory,
progress_file = "DM-predict.txt")## Starting predictions, running for 4000 iterations
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Bug - ESS for X shrunk to current position
## Prediction iteration 100
## Prediction iteration 200
## Prediction iteration 300
## Prediction iteration 400
## Prediction iteration 500
## Prediction iteration 600
## Prediction iteration 700
## Prediction iteration 800
## Prediction iteration 900
## Prediction iteration 1000
## Prediction iteration 1100
## Prediction iteration 1200
## Prediction iteration 1300
## Prediction iteration 1400
## Prediction iteration 1500
## Prediction iteration 1600
## Prediction iteration 1700
## Prediction iteration 1800
## Prediction iteration 1900
## Prediction iteration 2000
## Prediction iteration 2100
## Prediction iteration 2200
## Prediction iteration 2300
## Prediction iteration 2400
## Prediction iteration 2500
## Prediction iteration 2600
## Prediction iteration 2700
## Prediction iteration 2800
## Prediction iteration 2900
## Prediction iteration 3000
## Prediction iteration 3100
## Prediction iteration 3200
## Prediction iteration 3300
## Prediction iteration 3400
## Prediction iteration 3500
## Prediction iteration 3600
## Prediction iteration 3700
## Prediction iteration 3800
## Prediction iteration 3900
## Prediction iteration 4000n_samples <- length(pred$X[, 1])
N_pred <- length(pred$X[1, ])
X_ci <- apply(pred$X, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ]
caribou.df <- data.frame(
Covariate=c(X_ci),
Observation=factor(rep((1:N_pred), each=n_samples*0.95))
)
ggplot(caribou.df, aes(Observation, Covariate)) +
geom_violin(position="identity") +
scale_x_discrete(breaks=seq(5, 310, 5)) +
labs(x="Observation", y="Unobserved WTD")
y_prop <- y
y_recon_prop <- y_recon
for (i in 1:N) {
y_prop[i, ] <- y_prop[i, ] / sum(y_prop[i, ])
}
for (i in 1:N_pred) {
y_recon_prop[i, ] <- y_recon_prop[i, ] / sum(y_recon_prop[i, ])
}
## WA reconstruction - subset to deal with all zero occurrence species
zeros_idx <- which(colSums(y_prop) == 0)
if (length(zeros_idx) > 0) {
modWA <- rioja::WA(y_prop[, - zeros_idx], X)
predWA <- predict(modWA, y_recon_prop[, - zeros_idx], sse=TRUE, nboot=1000)
} else {
## no data to subset
modWA <- rioja::WA(y_prop, X)
predWA <- predict(modWA, y_recon_prop, sse=TRUE, nboot=1000)
}## Bootstrapping for SSE:
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======= | 12%
|
|======== | 12%
|
|======== | 13%
|
|========= | 13%
|
|========= | 14%
|
|========= | 15%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 19%
|
|============= | 20%
|
|============= | 21%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 28%
|
|=================== | 29%
|
|=================== | 30%
|
|==================== | 30%
|
|==================== | 31%
|
|==================== | 32%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|====================== | 35%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 36%
|
|======================== | 37%
|
|======================== | 38%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 45%
|
|============================== | 46%
|
|============================== | 47%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================= | 52%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 53%
|
|=================================== | 54%
|
|=================================== | 55%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================= | 61%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 62%
|
|========================================= | 63%
|
|========================================= | 64%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 70%
|
|============================================== | 71%
|
|============================================== | 72%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================ | 75%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|================================================== | 78%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 79%
|
|==================================================== | 80%
|
|==================================================== | 81%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 85%
|
|======================================================== | 86%
|
|======================================================== | 87%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|=========================================================== | 92%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================= | 95%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 96%
|
|=============================================================== | 97%
|
|=============================================================== | 98%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%pred_mu_WA <- predWA$fit[, 1]
pred_sd_WA <- sqrt(predWA$v1.boot[, 1])
## MLRC reconstruction - subset to deal with all zero occurrence species
zeros_idx <- which(colSums(y_prop) == 0)
if (length(zeros_idx) > 0) {
modMLRC <- rioja::MLRC(y_prop[, - zeros_idx], X)
predMLRC <- predict(modMLRC, y_recon_prop[, - zeros_idx],
sse=TRUE, nboot=1000)
} else {
modMLRC <- rioja::MLRC(y_prop, X)
predMLRC <- predict(modMLRC, y_recon_prop, sse=TRUE, nboot=1000)
}## Bootstrapping for SSE:
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======= | 12%
|
|======== | 12%
|
|======== | 13%
|
|========= | 13%
|
|========= | 14%
|
|========= | 15%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 19%
|
|============= | 20%
|
|============= | 21%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 28%
|
|=================== | 29%
|
|=================== | 30%
|
|==================== | 30%
|
|==================== | 31%
|
|==================== | 32%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|====================== | 35%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 36%
|
|======================== | 37%
|
|======================== | 38%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 45%
|
|============================== | 46%
|
|============================== | 47%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================= | 52%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 53%
|
|=================================== | 54%
|
|=================================== | 55%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================= | 61%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 62%
|
|========================================= | 63%
|
|========================================= | 64%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 70%
|
|============================================== | 71%
|
|============================================== | 72%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================ | 75%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|================================================== | 78%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 79%
|
|==================================================== | 80%
|
|==================================================== | 81%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 85%
|
|======================================================== | 86%
|
|======================================================== | 87%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|=========================================================== | 92%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================= | 95%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 96%
|
|=============================================================== | 97%
|
|=============================================================== | 98%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%pred_mu_MLRC <- predMLRC$fit[, 1]
pred_sd_MLRC <- sqrt(predMLRC$v1.boot[, 1])
## Modern analogue technique
modMAT <- mat(y, X)
predMAT <- predict(modMAT, y_recon, bootstrap=TRUE, n.boot=1000)
pred_mu_MAT <- apply(predMAT$predictions$model$predicted, 2, mean)
pred_sd_MAT <- apply(predMAT$predictions$model$predicted, 2, sd)
## GJAM model fit
idx_hold <- (N+1):(N+N_pred)
Xdf <- data.frame(x=c(X, rep(NA, N_pred)))
Xdf$x[idx_hold] <- NA
ydf <- data.frame(as.matrix(rbind(y, y_recon)))
colnames(ydf) <- paste("y", 1:dim(y)[2], sep="")
ml <- list(ng = 5000, burnin = 500, typeNames = rep("CC", dim(y)[2]),
PREDICTX=TRUE)
## fit second order polynomial model
# out <- gjam(~ x + I(x^2), Xdf, ydf, ml)
#
# pred_mu_GJAM <- out$prediction$xpredMu[idx_hold, 2] #inverse prediction of x
# pred_sd_GJAM <- out$prediction$xpredSd[idx_hold, 2]
## Random Forest
library(randomForest)## randomForest 4.6-12## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margintrain <- data.frame(y=X, as.matrix(y))
test <- data.frame(y_recon)
n_samples <- length(pred$X[, 1])
rf <- randomForest(y ~ ., data = train, ntree=n_samples)
preds_rf <- predict(rf, test, predict.all=TRUE)$individualn_samples <- length(pred$X[, 1])
N_pred <- length(pred$X[1, ])
preds_WA <- matrix(0, n_samples, N_pred)
preds_MLRC <- matrix(0, n_samples, N_pred)
preds_MAT <- matrix(0, n_samples, N_pred)
# preds_GJAM <- matrix(0, n_samples, d)
for (i in 1:N_pred) {
preds_WA[, i] <- rnorm(n_samples, pred_mu_WA[i], pred_sd_WA[i])
preds_MLRC[, i] <- rnorm(n_samples, pred_mu_MLRC[i], pred_sd_MLRC[i])
preds_MAT[, i] <- rnorm(n_samples, pred_mu_MAT[i], pred_sd_MAT[i])
# preds_GJAM[, j] <- rnorm(n_samples, pred_mu_GJAM, pred_sd_GJAM)
}
X_ci <- rbind(apply(pred$X, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_WA, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_MLRC, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_MAT, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ],
apply(preds_rf, 1, sort)[(0.025*n_samples+1):(0.975*n_samples), ])
# apply(preds_GJAM, 2, sort)[(0.025*n_samples+1):(0.975*n_samples), ])
caribou.df <- data.frame(
Covariate=c(X_ci),
Observation=factor(rep((1:N_pred), each=5*n_samples*0.95)),
Model=rep(c("MVGP", "WA", "MLRC", "MAT", "RF"), each=n_samples*0.95)
)
ggplot(caribou.df, aes(Observation, Covariate, color=Model)) +
geom_violin(position="dodge") +
scale_x_discrete(breaks=seq(5, 310, 10)) +
labs(x="Observation", y="Unobserved Water Depth") +
facet_wrap(~ Model, ncol=2) +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model))
ggsave("caribou_recons.png")## Saving 7 x 5 in imageggplot(caribou.df, aes(Observation, Covariate, color=Model)) +
#geom_violin(position="dodge") +
scale_x_discrete(breaks=seq(310, 0, -10)) +
labs(x="Observation", y="Unobserved Water Depth") +
facet_wrap(~ Model, ncol=2) +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
scale_y_reverse()
ggplot(caribou.df, aes(Observation, Covariate, color=Model)) +
labs(x="Observation", y="Unobserved Water Depth") +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
stat_summary(fun.ymin = function(z) {quantile(z, 0.025)},
fun.ymax = function(z) {quantile(z, 0.975)},
geom = "ribbon", aes(Observation, Covariate,
group=Model, color=Model, fill=Model), alpha=0.25) +
ylim(-10, 60)## Warning: Removed 15810 rows containing non-finite values (stat_summary).
## Warning: Removed 15810 rows containing non-finite values (stat_summary).
ggplot(subset(caribou.df, Model %in% c("MVGP", "WA")), aes(Observation, Covariate, color=Model)) +
labs(x="Observation", y="Unobserved Water Depth") +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
# stat_summary(fun.ymin = function(z) {quantile(z, 0.025)},
# fun.ymax = function(z) {quantile(z, 0.975)},
# geom = "ribbon", aes(Observation, Covariate,
# group=Model, color=Model, fill=Model), alpha=0.25) +
ylim(-10, 60)## Warning: Removed 230 rows containing non-finite values (stat_summary).
ggplot(caribou.df, aes(Observation, Covariate, color=Model)) +
labs(x="Observation", y="Unobserved Water Depth") +
stat_summary(fun.y = mean, geom = "line", aes(Observation, Covariate, group=Model)) +
stat_summary(fun.ymin = function(z) {quantile(z, 0.025)},
fun.ymax = function(z) {quantile(z, 0.975)},
geom = "ribbon", aes(Observation, Covariate,
group=Model, color=Model, fill=Model), alpha=0.25) +
facet_wrap(~Model, ncol=2) +
ylim(-10, 60)## Warning: Removed 15810 rows containing non-finite values (stat_summary).
## Warning: Removed 15810 rows containing non-finite values (stat_summary).